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Abstract 

We propose a new efficient scheme for the quantum Monte Carlo study of 
quantum critical phenomena in quantum spin systems. Rieger and Young's 
Trotter-number-dependent finite-size scaling in quantum spin systems and Ito 
et aVs evaluation of the transition point with the non-equilibrium relaxation in 
classical spin systems are combined and generalized. That is, only one Trotter 
number and one inverse temperature proportional to system sizes are taken 
for each system size, and the transition point of the transformed classical spin 
model is estimated as the point at which the order parameter shows power-law 
decay. The present scheme is confirmed by the determination of the critical 
phenomenon of the one-dimensional 5 = 1/2 asymmetric XY model in the 
ground state. The estimates of the transition point and the critical exponents 
(3, 7 and v are in good agreement with the exact solutions. The dynamical 
critical exponent is also estimated as A = 2.14 ±0.06, which is consistent with 
that of the two-dimensional Ising model. 

1 Introduction 

Quantum critical phenomena are one of the most interesting topics in low- dimensional 
quantum spin systems, and various numerical methods have been proposed. In one- 
dimensional systems, the exact-diagonalization method gives the most reliable result 
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when the chain is single and the size of spins is not so large. In multi-leg ladder 
systems or large-spin systems, reachable system sizes with this method are not sat- 
isfactory and other approximate methods should be used. Historically, the quantum 
Monte Carlo (QMC) method]!], [J and related methods (the quantum transfer-matrix 
method p, |3], f|, the Monte Carlo power method 0, |[ etc.) have been used, and nowa- 
days the density- matrix renormalization-group (DMRG) method [|, ^, || is considered 
to be most powerful because quite large (effectively infinite) systems can be calcu- 
lated with small errors and relatively small required CPU demand. The DMRG 
method is applicable even to two-dimensional systems with large enough energy gap, 
but the limitation of computer memories makes the estimation of critical phenomena 
difficult at present (any systems are gapless at the critical point in the thermody- 
namic limit). Then, the QMC method is still most useful to study quantum critical 
phenomena in non-frustrated two-dimensional quantum spin systems. 

In QMC simulations of ground states in quantum spin systems, the most serious 
problem is large sampling errors. The QMC method is essentially applicable only to 
finite-temperature systems, and "ground states" mean nothing but "calculations at 
low enough temperatures" , which might include ambiguity. Moreover, simulations at 
low temperatures result in large correlation time, which is the origin of large sampling 
errors. Furthermore, in the conventional world-line QMC method, physical quantities 
are evaluated after the extrapolation of the data for various Trotter numbers, which 
makes statistical errors larger. These shortcomings of the world-line QMC algorithm 
have been overcome with the stochastic series expansion method|| 10, 11], the quan- 
tum cluster algorithm [12, [13| or Prokof 'ev's algorithm 0, |ll| proposed recently. The 
Trotter extrapolation is not required in these methods, and fast relaxation of them 
makes simulations at very low temperatures possible. As long as finite systems are 
considered, energy gaps are finite even at the transition point, and simulations at 
the temperatures lower enough than energy gaps give "ground-state" properties with 
convergent correlation time. Although such a brute-force strategy is required for the 
evaluation of actual physical quantities at the ground state, more efficient scheme is 
expected to study quantum critical phenomena, in which universal behavior exists. 

In the present paper, we propose a new efficient scheme to study quantum critical 
phenomena in quantum spin systems with the QMC method. This scheme can be 
regarded as a combination and a generalization of two independent techniques of 
Monte Carlo calculations originally proposed in Ising spin systems. That is, Rieger 
and Young's Trotter-number-dependent finite-size scaling[|l6|] introduced in the QMC 
study of the quantum critical phenomenon of the Ising spin glass model in a trans- 
verse field is coupled with Ito et. a/.'s determination of the transition point |I7|] based 
on the non-equilibrium relaxation of the order parameter. The latter method was 
introduced in the classical Ising spin glass model, and the former method is necessary 
for the generalization to quantum spin systems. By use of the present scheme, only 
one Trotter number has to be considered in each system size for the determination 
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of the universality class. In short, the combination of these two methods is essential 
for efficient numerical calculations in quantum spin systems. 

In order to verify the applicability of the present scheme, we estimate the tran- 
sition point and critical exponents of the one-dimensional 5 = 1/2 asymmetric XY 
model. This model is mapped on a two-dimensional Ising spin model with conserva- 
tion laws (or infinite multi-body interactions in the Ising-spin representation) by the 



Suzuki- Trotter decomposition!], and the universality class of this model Pq, [19| , |20 
is different from that of the two-dimensional Ising model. In section 2, the present 
scheme is explained together with a brief review of the world-line QMC method. 
In section 3, the quantum critical phenomenon of the one-dimensional 5 = 1/2 
asymmetric XY model is analyzed. The transition point and the dynamical critical 
exponent are estimated with very large clusters (linear sizes are L = 1024, 2048 and 
4096) and various critical exponents are estimated with the clusters up to L = 128. 
The above descriptions are summarized in section 4. Part of the present results have 



already been briefly reported elsewhere [21 



2 Numerical methods 

In the world-line quantum Monte Carlo method]!], [|, the partition function of quan- 
tum spin systems is transformed to a classical one by the Suzuki- Trotter formula, 

r fa \ la \1 M 




exp — — rti] ex P 




where (3 stands for the inverse temperature {(3 = 1/T), and the partial Hamiltonians 
Tii and H.2 in one-dimensional systems are shown in figure [|, which is known to be 
the checkerboard decomposition. Then, each 4-spin unit cell is independent, and 
Monte Carlo simulations are possible. Conventional local flips and the x-direction 
global flip are introduced. When the models have conservation laws, the Trotter- 
direction global flip should be considered for the evaluation of the thermal averages 
in the whole phase space. As is well known, the tangent-direction global flip can be 



safely omitted, and this algorithm is vectorized |22[ . 

In the conventional world-line QMC method, quantities in the ground state are 
evaluated as follows: 

1. The temperature is fixed at a "low enough" one, T (L), which is much smaller 
than the energy gap of finite clusters, AE(L), in the subspace in which the 
ground state is included. The Trotter-direction global flip is not introduced in 
this case. 

2. QMC simulations are performed for some finite Trotter numbers M, and quan- 
tities in the M — > oo limit are estimated with the following scaling form, 
Q(M, L) = Q(oo, L) + const. /M 2 + ■■■. 
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3. Quantities in the thermodynamic limit are evaluated from the finite-size scaling 
of {Q(oo, L)} for various system sizes. 

In the above analysis, the T — > limit is skipped by use of the finite energy gap in 
finite systems, but double extrapolations with respect to M and L are still neces- 
sary. Calculations for some Trotter numbers are required for each system size, and 
simulations at a "low enough" temperature consume much CPU demand because of 
long correlation time. 

Although the above analysis is required for the evaluation of a certain physical 
quantity in the ground state, such a complicated procedure is not necessary for the 
investigation of quantum critical phenomena. The universality class of the original 
D-dimensional quantum system in the ground state is expected to be equivalent 
to that of the transformed (D + l)-dimensional classical system, and Rieger and 
Young pointed out that the finite-size scaling holds when the ratios M/L z and 
(3/M are fixed. The exponent z stands for the ratio of the correlation length in the 
real direction (£) and that in the Trotter direction (£ T ), namely £ z oc £ T in the vicinity 
of the transition point. This exponent is given by z = 1 in non-random systems (the 
Lorentz invariance). In this scheme, the temperature is scaled in accordance with 
the system size, and the ambiguity of the simulated temperature is automatically 
removed. Moreover, this temperature can be higher than the standard "low enough" 
temperature. 

The transition point has been estimated from physical quantities in equilibrium 
until present. When it is determined by the finite-size scaling analysis, the critical 
exponent should also be evaluated simultaneously, and the accuracy of estimates 
becomes poor. When it is estimated from scale-invariant quantities, the above dis- 
advantage is excluded. However, statistical errors of such complicated quantities are 
larger than those of ordinary quantities, and the final estimates are not expected 
to be so accurate again. Then, in the present paper, we utilize a new method for 
the determination of the transition point recently proposed by Ito et. al.\\L7\. They 
showed that the transition point can be estimated quite accurately with the non- 
equilibrium relaxation of the order parameter. That is, simulations are started from 
a perfectly-ordered state, and the order parameter is measured in each Monte Carlo 
step (MCS). The order parameter decays exponentially in the disordered phase, re- 
mains non-vanishing in the ordered phase, and power-law decay is observed only 
at the transition point (figure This calculation uses only several hundreds or a 
few thousands of MCS in each sample, and very large systems can be easily treated 
without equilibration and Monte Carlo sampling. 

Since the Trotter extrapolation is justified only in equilibrium quantities, this 
method is not consistent with the conventional world-line QMC method. Thus, 
the combination with Rieger and Young's scheme is indispensable for applications 
to quantum spin systems. Since the transition point is not a universal value, the 
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estimate obtained from transformed classical spin models with fixed (3/M is generally 
different from the transition point of the original quantum spin model. However, the 
transition point of the transformed model should be used in the finite-size-scaling 
analysis of a series of systems with fixed /3/M, and the "true" transition point of 
the original model is estimated from some simulations with different values of /3/M. 
Details of this analysis will be explained elsewhere ||23|| . 

The present scheme is so general that it can be coupled with any kinds of Monte 
Carlo dynamics in principle. Although values of the dynamical critical exponent 
depend on the dynamics (details will be explained in the next section), the power- 
law-decay behavior at the transition point is common as long as second-order phase 
transitions are considered. The essential point of the usage of Rieger and Young's 
scheme in the present framework is the idea that simulated temperatures are scaled 
in accordance with variance of system sizes. Thus, Trotter-number-free methods such 
as the continuous-time loop algorithmfl3| can be utilized. On that occasion, j3/L z 
should be fixed instead of M/L z and j3/M. On the other hand, the accuracy of the 
estimate of the transition point depends on the dynamics. If a very fast dynamics 
such as the loop algorithm is utilized, the decay process is completed within first 
several MCS and precise evaluation of the transition point would be difficult because 
of large fluctuations of the order parameter and the discontinuity of MCS. Of course, 
this statement is also model dependent: even the loop algorithm might be useful in 
systems with principal slow dynamics such as quantum spin glass models. We can 
at least tell that the conventional "slow" dynamics is suitable for regular systems as 
will be shown in the next section. 



3 Numerical calculations of critical phenomenon 

In order to verify the applicability of the new scheme explained in the previous 
section, we analyze the quantum critical phenomenon of the one-dimensional S = 1/2 
asymmetric XY model described by the following Hamiltonian, 

w = -E S ? S t+i " 9 E s?s? +1 -H X J2 s? , (2) 

i=l i=l i=l 

with the periodic boundary condition = Si' y ). The in-plane field H x is intro- 

duced only for the definition of the zero-field susceptibility, and it is not applied in 
actual simulations. We only consider even-spin systems which have a singlet ground 
state. This is the simplest quantum spin model with conservation laws, and this 
model has already been solved exact ly [|T8|, |19| , [2~0 |. There are no matrix elements 
between the two subspaces with J2i Sf = odd and J2i Sf = even. The apparent tran- 
sition point g c = 1 is common in the original model (0) and the transformed classical 
spin models. Simulations are based on the up-down basis along the x-direction, 
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and the Trotter number and simulated temperature are fixed as 2M / L = 1 and 



{3/M = 1, following Rieger and Young's parameterization [jlB] . Under these condi- 



tions, the original L-spin system is transformed to a classical L x L-spin system with 
all the interactions unity at g — 1, and good scaling behavior is expected. The phase 
space is limited to the one with J2i Sf = even in which the ground state is included. 
Since finite (not "low enough") temperatures are introduced in the present scheme, 
the above restriction of the phase space cannot be justified a priori. Complementary 
simulations without this restriction are discussed in Appendix. 



3.1 Transition point g c 

First, the transition point is estimated with the non-equilibrium relaxation of the 
order parameter. When the x-component of the magnetization, 

m*(L) = jJ2S? , (3) 

i=l 

is regarded as the order parameter, the regions g < 1 and g > 1 correspond to 
the ordered and disordered phases, respectively. Simulations are started from the 
all-up state, and the initial 3 x 10 3 Monte Carlo steps (MCS) are measured for the 
L = 1024, 2048 and 4096 systems. In simulations for these system sizes, data with 
20, 5 and 1 independent sequences of random numbers are averaged, respectively. 
The order parameter of the L = 4096 system is plotted versus MCS in a logarithmic 
scale in figure f|(a), and this behavior is the same as displayed in figure |j. That is, 
the data at g = 1.000 are almost on a straight line, and the data at g = 0.996 and 
g = 1.004 are clearly upward- and downward-bending, respectively. Therefore, the 
transition point is estimated as g c = 1.000 ± 0.002 within the present simulations, 
as expected. The data for the L = 1024 and 2048 systems at g — 1.000 are given 
in figure |](b) together with those for the L = 4096 system, and this figure means 
that there exists no explicit systematic size dependence in such large systems. Since 
values of I x 2M spins are averaged in these simulations, the data show relatively 
smooth behavior in spite of small number of averaged random-number sequences. 

Next, we perform similar simulations in smaller systems. In equilibrium calcula- 
tions, the "approximate" transition point ga, p (L) such as the maximum of the specific 
heat shows the following systematic deviation from the transition point g c in finite 
systems, 

g ap {L) =g c + const, x L 1/v + ■■■ , (4) 

and the critical exponent of the correlation length, z/, is estimated. The relaxation of 
the order parameter for L = 32 and 64 (the numbers of averaged samples are 16384 
and 4096, respectively) is displayed in figures |](a) an d H(b), respectively. These 
figures show that the order parameter decays exponentially after a certain MCS 
even for g < 1. This property is quite natural because the order parameter vanishes 
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regardless of the value of g in finite systems in equilibrium. These figures also reveal 
that the crossover MCS from "quasi-infinite" properties to finite-size ones depends on 
system sizes, namely ~ 40 MCS for L = 32, and ~ 180 MCS for L = 64. These results 
suggest that the above behavior of L > 1000 systems is nothing but transient, but 
that the present 3 x 10 3 MCS are small enough to be regarded as "quasi-infinite" time 
scale. Moreover, the transition point estimated from these figures is also g c ~ 1.00, 
and the size dependence of the estimates seems quite small. This behavior is not 
unphysical because the "quasi-infinite" property means that finite-size effects are not 
apparent in this time scale. Thus, the scaling form (Q) is not useful for the estimation 
of the critical exponent v within the present scheme. 

As explained in the previous section, the "transition point" estimated with the 
present scheme generally depends on the value of /3/M. On the other hand, it is 
independent of the value of M/L, because "quasi-infinite" properties are observed 
in the present scheme and the difference of the shape of transformed clusters from 
squares is expected to result only in the difference the crossover MCS as shown above. 
Then, in order to confirm this argument, we simulate systems with (a) 4M/L = 1 
(L = 4096, M = 1024, 2 samples) and (b) M/L = 1 (L = 2048, M = 2048, 2 
samples). These results are plotted versus MCS in a logarithmic scale in figures |](a) 
and Mb), respectively. These figures also give g c = 1.000 ± 0.002, as expected. 



3.2 Dynamical critical exponent A 



According to the dynamical finite-size scaling theory f24], g5| , the non-equilibrium 
relaxation of the order parameter at the transition point is scaled [^6| as 



ro*(t)~r\ A = A (5) 
Av 

where t stands for MCS, the critical exponents (3 and v will be defined in the next 
subsection ((|i~4|) and ([TBI) ), and the dynamical critical exponent A is defined as 

(Smsm) - (Sf (0)> 2 ~ expH/r) , (6) 
T~\g-g c \~ A for g -> g c ■ (7) 

Although the transition point can be estimated accurately enough with the data 
displayed in figure |](a), these data are too divergent to evaluate the exponent A with 



the local-exponent method originally proposed by Stauffer|26|. Then, we naively 
fit the data shown in figure |4](a) with the scaling form for various time scales 
{t = ti n i ~ tfin, At = tfi n — i ini = 50, 100 and 200). Estimates of the exponent A 
are plotted versus 1 / t m id i n figure £|, with t mid = (t ini + tfi n )/2. In addition to the 
data given in figure |](a), results of two more independent simulations are averaged 
in order to obtain reliable estimates. In each At, A almost remains constant in 
the early stage, and it begins to fluctuate as t mid becomes larger. As for the final 
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estimate A est , we require that the central value should be close to the early-stage one, 
and that initial up-down fluctuations should be included within the error bar. These 
conditions are satisfied when we take A cs t = 0.1167±0.0033 (see figure f|). Combining 



this exponent with the exact solution, (3ju — 0.25 JT^, [19], p0| , the dynamical critical 
exponent of the present model is estimated as 

A = 2.14 ±0.06, (8) 

which is consistent with that of the two-dimensional Ising model, A = 2.165 ± 



0.010 p8f. This result might suggest that the dynamical critical exponent is com- 
mon in the free-fermion models, even though the universality classes are different. 
Actually, the critical exponent v defined by 

(^)-(^) 2 ~exp(-|r,-r,|/0 , (9) 

€~\g-9c\~ v for 9 -> 9c , (10) 

is unity in both the present model and the two-dimensional Ising model. When the 
definition of A, @ and (0), is compared with that of u, @ and fllOD , the exponent 
A can be regarded as the MCS version of v. Since these two models have the 
common u, it is not strange that they might also have the common A. Of course, 
the statistical error of the present estimate is not small enough, and further studies 
in this direction (averaging much sequences of random numbers or simulating larger 
systems) are required. 

3.3 Critical exponents (3/v and 7/1/ 

Next, we estimate the critical exponents which can be obtained from equilibrium 
simulations at g = g c = 1. We evaluate the squared magnetization defined by 

m ' to ,L) S (<m-(L))») = ^KfgS^) 2 ) , (11) 
and the zero-field susceptibility for the in-plane field defined by 

L 2M 



x{9,L)=jt-(m'(L)) 







HX=0 4M*L 



L 2M \"\ 
i=l n =l J I 



(12) 



Here n stands for the suffix in the Trotter direction, and these quantities show the 
following scaling behavior at the transition point, 

m 2 (g c ,L)~L-W , X (g c , L) ~ U' v , (13) 

where the critical exponents /3, 7 and v are defined as 

m(g, 00) ~ (g c - gf for g -> g c , (14) 

x(g,oo) ~ Ig-gd^ for g -> g c , (15) 

f(flr,oo) ~ l^-^cr" for 9 9c ■ (16) 



8 



Table 1: Monte Carlo steps for equilibration (E-MCS) and for measurement (M- 
MCS) are listed for various system sizes. 



L 


E-MCS 


M-MCS 


16 


0.2 x 10 6 


0.2 x 10 7 


24 


0.4 x 10 6 


0.4 x 10 7 


32 


0.8 x 10 6 


0.8 x 10 7 


48 


1.2 x 10 6 


1.2 x 10 7 


64 


0.5 x 10 6 


0.5 x 10 7 


80 


0.8 x 10 6 


0.8 x 10 7 


96 


1.0 x 10 6 


1.0 x 10 7 


112 


1.2 x 10 6 


1.2 x 10 7 


128 


2.0 x 10 6 


2.0 x 10 7 



Simulations are performed for the L = 16, 24, 32, 48, 64, 80, 96, 112 and 128 systems, 
and the MCS for each system size are listed in table [l]. Results of these calculations 
are given in table |2| and figures ^(a) and [|(b). The critical exponents are estimated 
by the least-squares fitting of the data for 64 < L < 128 with the scaling forms ([13]) 

as 

p/v = 0.241 ± 0.003 , 7 /z/ = 1.499 ± 0.005 . (17) 

On the other hand, when we utilize the following scaling forms in which the next- 
order correction terms are considered, 

m 2 (g c ,L) = a.L- 2 ^ + ^L- 2 ^- 1 + ■ ■ ■ , (18) 
X {g c ,L) = bl L^ + b^- 1 + ■ ■ ■ , (19) 

we have 

P/u = 0.245 ± 0.006 , 7/z/ = 1.49 ± 0.02 , (20) 
from the least-squares fitting of the data for 24 < L < 128. The estimates (O) and 



(GO) are both consistent with the exact values, (3/v = 0.25 and 7/z/ = 1.5 P^, |T9|, |20| 



3.4 Critical exponent v 



The critical exponents (3/v and 7/z/ are related[29] to one another through the ex- 
ponent z as 

2/3/t/ + 7/z/ = d + z , (21) 

where d denotes the spatial dimension. Then, the critical exponent v should be 
estimated independently for the complete determination of the universality class. 



9 



Table 2: Ground-state averages of the squared magnetization and the susceptibility 
at the transition point for various system sizes. 



L 


m 2 


16 


0.11094 ±0.00031 


24 


0.09379 ± 0.00035 


32 


0.08279 ±0.00037 


48 


0.06868 ± 0.00023 


64 


0.06052 ± 0.00029 


80 


0.05437 ±0.00021 


96 


0.04982 ± 0.00026 


112 


0.04624 ± 0.00026 


128 


0.04336 ± 0.00022 



X 

1.2242 x 10 1 ±3.0 x 10~ 2 
2.2785 x 10^6.5 x 10" 2 
3.531 x 10 1 ± 1.2 x 10- 1 
6.589 x 10 1 ±3.0 x 10" 1 
1.0299 x 10 2 ±6.3 x 10" 1 
1.4385 x 10 2 ± 7.2 x lO^ 1 
1.890 xl0 2 ± 1.3x10° 
2.381 xl0 2 ± 1.7x10° 
2.912 xl0 2 ± 1.9x10° 



This calculation requires several simulations out of the transition point, and they 
consume much more CPU time than the simulations in the previous subsection. 
Since the direct evaluation of the correlation length is not so accurate, we instead 
estimate the critical exponent v on the basis of the following scaling functions, 

m 2 (g,L) ~ L-^f(L/ag))^L- 2 ^F(L^\g-g c \), (22) 
x (g,L) ~ L->'»g(L/t(jg)) ~ L^G(L^\g - g c \) . (23) 

These formulas mean that the derivatives with respect to g behave as 
d 



-m 2 (g,L) 

dg 



~L^»», ± X (g,L) 
dq 

a=9c ^ 



9=9c 



Practically, such derivatives are replaced by the evaluation of the slopes of the phys- 
ical quantities just below the transition point. For example, the susceptibility of 
the L = 128 system is displayed in figure |j. The exponents appeared in ([24]) are 
estimated as 

(1 -2p)/v = 0.525 ±0.007 , (1 + 7)/// = 2.506 ± 0.006 , (25) 

by the least-squares fitting of the data for 80 < L < 128 (see table || and figures |(a) 
and H(b)) with (p4|) . From the exponents (fP7|) and (|25|), we have 

l/i/ = 1.007 ±0.008 or v = 0.993 ± 0.008 , (26) 

assuming the statistical independence of the exponents (|17|) and (|25|) . This estimate 
is consistent with the exact value, v = l|pO|l. Note that the complete coincidence of 
the estimates of v obtained from two different quantities would be by chance. 
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Table 3: Slopes of the squared magnetization and the susceptibility just below the 
transition point for various system sizes. 



L 


m 2 


X 


64 


0.5942 ± 0.0039 


1301 ± 8 


80 


0.6804 ± 0.0062 


2320 ± 21 


96 


0.7590 ± 0.0086 


3720 ± 42 


112 


0.8191 ±0.0086 


5457 ±57 


128 


0.8714 ±0.0074 


7567 ± 64 



4 Summary and future problems 

In the present paper, we have proposed a new efficient finite-size scaling scheme 
to study quantum critical phenomena in quantum spin systems on the basis of the 
quantum Monte Carlo method. Instead of the conventional Trotter extrapolation 
at low enough temperatures, we concentrate on transformed Ising spin systems with 
2M/L = 1 and /3/M = 1. The transition point of the transformed model is estimated 
with the non-equilibrium relaxation of the order parameter, and various critical ex- 
ponents are evaluated from the standard finite-size scaling of transformed classical 
spin systems. That is, the critical exponents divided by the critical exponent of the 
correlation length, z/, are estimated from the size dependence of physical quantities 
at the transition point. The exponent v itself is evaluated from the size dependence 
of the slopes of physical quantities just below the transition point. Then, the univer- 
sality class of the models are completely determined within the present framework. 
In order to verify the applicability of the present scheme, we have estimated the tran- 
sition point and the critical exponents f3, 7 and v of the one-dimensional S — 1/2 
asymmetric XY model, and all the critical exponents have been evaluated with high 
accuracy in comparison with the exact solutions. The dynamical critical exponent 
has also been estimated as A = 2.14 ± 0.06, which is consistent with that of the 
two-dimensional Ising model. The present estimation of the transition point is more 
efficient than any other methods based on equilibrium quantities. The present finite- 
size-scaling scheme is so general that it can be combined with improved algorithms 
such as the quantum cluster algorithm. 

The present simulations have been performed in the restricted phase space with 
J2i Sf = even in which the ground state is included. Complementary simulations 
without this restriction (simulations including the Trotter-direction global flip) have 
also been performed, and we have found that theses simulations require much CPU 
demand (about 1.3 times larger in the present model) and show larger finite-size 
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corrections than those with the restriction of the phase space. These results mean 
that simulations including the Trotter-direction global flip are useless. We emphasize 
that the conventional "slow" single-spin-flip-type dynamics is practically useful for 
the precise estimation of the transition point in the present scheme. On the other 
hand, the introduction of the quantum cluster algorithm in simulations of equilib- 
rium quantities is an important future task for us. Applications to quantum critical 
phenomena of the two-dimensional 5=1/2 dimerized Heisenberg models |23j and 
the two-dimensional 5* = 1/2 random Heisenberg models^] will be reported in the 
near future. 
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Appendix 

Here we show the results of simulations including the Trotter-direction global flip. 
Since details of the analysis have already been explained in section 2 and 3, only the 
results for the transition point and the critical exponents (3/v and j/u are briefly 
summarized. 

The transition point is estimated with the L = 4096 system for one sequence of 
random numbers. The decay of the order parameter from the all-up state is measured 
during 3 x 10 3 MCS as shown in figure |], and we have g c = 1.000 ± 0.002. Although 
this estimate is the same as obtained in section 3, fluctuations of the order parameter 
become larger than those measured in section 3. This result means that the present 
scheme also works well when the Trotter-direction global flip is included, but the 
expansion of the phase space results in the increase of fluctuations, which is not 
desirable for accurate estimation of critical phenomena. 

The squared magnetization and the zero-field susceptibility are simulated at the 
transition point g = g c = 1 for the same system sizes and the same MCS as listed 
in table [p. The size dependence of these quantities are displayed in figures f|(a) and 
[|(b), and the critical exponents are estimated with the scaling forms ( |TB| ) as 

p/v = 0.219 ±0.007 , 7/z/ = 1.546 ± 0.009 , (27) 
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from the 64 < L < 128 clusters, and with (|i~8|) and ([19]), we have 

/3/z/ = 0.22 ± 0.03 , 7/1/ = 1.56 ±0.05 , (28) 

from the 32 < L < 128 clusters. These estimates are not inconsistent with the exact 
solutions, (3/v = 0.25 and 7/1/ = 1.5, but the accuracy is not so good as the results 
obtained in section 3 due to larger finite-size corrections. 

The above two numerical calculations with respect to non-equilibrium relaxation 
and equilibrium quantities are based on the same algorithm, and they consume larger 
CPU time than those in section 3 because of the extra Trotter-direction global flip. 
They spend about 1.3 times larger CPU demand in the present model. All the results 
suggest that simulations including the Trotter-direction global flip are inferior to 
those without that kind of global flip. 
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Figure captions 

Figure 1: The checkerboard decomposition in one-dimensional systems. 

Figure 2: Schematic behavior of the non-equilibrium relaxation of the order pa- 
rameter from a perfectly-ordered state. Values of the order parameter are plotted 
versus Monte Carlo steps in a logarithmic scale (a) in the ordered phase, (b) at the 
transition point, and (c) in the disordered phase. 

Figure 3: The order parameter of the present model is plotted versus Monte Carlo 
steps in a logarithmic scale (a) for L = 4096 at g — 0.996 ~ 1.004 and (b) for 
L = 1024, 2048 and 4096 at g — 1.000. Solid lines are drawn as visual guides in the 
first figure. 

Figure 4: The order parameter of the present model is plotted versus Monte Carlo 
steps in a logarithmic scale for (a) L = 32 and (b) L = 64. The crossover from 
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"quasi-infinite" properties to finite-size properties is clearly seen. The crossover time 
scale is shown by the dashed line in each figure. 

Figure 5: The order parameter of the present model is plotted versus Monte Carlo 
steps in a logarithmic scale (a) for L = 4096 and M = 1024 at g = 0.996 ~ 1.004 
and (b) for L = 2048 and M = 2048 at g = 0.996 ~ 1.004. Solid lines are drawn as 
visual guides. 

Figure 6: Estimates of the exponent A for various ranges of time scales (t = t ini ~ t fin , 
At = t fin — tini = 50, 100 and 200) are plotted versus l/t m id, with t mid = (tmi + t fin )/2. 
The solid, dotted and bold-solid curves correspond to the results for At = 50, 100 
and 200, respectively. The final estimate A cs t = 0.1167 ± 0.0033 is displayed by 
straight lines. 

Figure 7: Ground-state averages of (a) the squared magnetization and (b) the suscep- 
tibility at the transition point are plotted versus system sizes in a logarithmic scale. 
The solid lines are drawn by the least-squares fitting of the data for 64 < L < 128 
with (11), and the dashed lines are obtained by the least-squares fitting of the data 
for 24 < L < 128 with (16) and (17). 

Figure 8: The susceptibility in the L = 128 cluster is plotted versus the control 
parameter g in the vicinity of the transition point. The slope of these data (the solid 
line) is evaluated by the least-squares fitting of the data drawn with full symbols. 

Figure 9: The slopes of (a) the squared magnetization and (b) the susceptibility are 
plotted versus system sizes in a logarithmic scale. The solid lines are drawn by the 
least-squares fitting of the data for 80 < L < 128 with (22). 

Figure 10: The order parameter of the present model including the Trotter-direction 
global flip is plotted versus Monte Carlo steps in a logarithmic scale for L = 4096. 
Solid lines are drawn as visual guides. 

Figure 11: Size dependence of (a) the squared magnetization and (b) the susceptibil- 
ity at the transition point in a logarithmic scale. The Trotter-direction global flip is 
included in these simulations. The solid lines are drawn by the least-squares fitting 
of the data for 64 < L < 128 with (11), and the dashed lines are obtained by the 
least-squares fitting of the data for 32 < L < 128 with (16) and (17). 



15 



16 




►real 



Fig.l: Y.Nonomura 



Fig. 2: Y.Nonomura 



(a) 
(b) 
(c) 

►log(MCS) 



Fig. 3(a): Y.Nonomura 




012345678 

log(MCS) 



Fig. 3(b): Y.Nonomura 



~ i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — r 



■0.8 

o L=1024 



1 .0 

1=2048 

E " 1 - 2 [ V ° L=4096 

O) 

O -1 .4 



1.6 - 



1 .8 



(b) 



012345678 

log(MCS) 



Fig. 4(a): Y.Nonomura 




□- 

j i i i i i i i i i i i i i i i i i_i i i i i i i i i_ 



1 2 3 4 5 

log(MCS) 



Fig. 4(b): Y.Nonomura 




Fig. 5(a): Y.Nonomura 




j i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i EG 



1 2 3 4 5 6 7 

log(MCS) 



Fig. 5(b): Y.Nonomura 




012345678 



log(MCS) 



Fig. 6: Y.Nonomura 



■0.105 



■0.110 



c< -0.1 15 - 



■0.120 



■0.125 - 



At= 50 
At = 100 
At=200 




mid 



Fig. 7(a): Y.Nonomura 




Fig. 7(b): Y.Nonomura 




Fig. 8: Y.Nonomura 




Fig. 9(b): Y.Nonomura 




Fig. 10: Y.Nonomura 




01 2345678 

log(MCS) 



Fig.11(a): Y.Nonomura 




Fig. 11(b): Y.Nonomura 




